
**************************************************************************************
**************************************************************************************
********************************** MAIN FIGURES **************************************
**************************************************************************************
**************************************************************************************

* NOTE: All file paths and data names that need to be inserted are identified with <>

clear
set more off
capture log close 

global dta = "<insert path to raw data>" 
global output= "<insert path to outputs>" 

log using "${output}\log\create_figures.log", replace

di "Time $S_DATE $S_TIME"

/* Programs that need to be installed:
rdplot 
regsave 
*/

**************************************************************************************

**************************************************************************************
************************************ Figure 2 ****************************************
**************************************************************************************
*** Mortality Rates around the Bend Points. (Produces Panels A, B & C.)

local sample main

**************************************************************************************
*** A: Lower bend point

*** Local macros
local bp lower 		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "d_year14"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`sample'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`sample'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-400(200)400"
local ylabel "3.4(0.3)4.4"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/Fig2A.eps", replace

**************************************************************************************
*** B: Family bend point

*** Local macros
local bp fm 		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`sample'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`sample'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "3.1(0.3)4.1"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/Fig2B.eps", replace

**************************************************************************************
*** C: Upper bend point

*** Local macros
local bp upper		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`sample'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`sample'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Upper bend point"
local yaxis "Mortality rate (%)"
local xaxis "Distance from Bend Point ($)"
local xlabel "-600(200)600"
local ylabel "4.4(0.2)5.1"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/Fig2C.eps", replace


**************************************************************************************
************************************ Figure 3 ****************************************
**************************************************************************************
*** Mortality Estimates using Varying Bandwidths.

**************************************************************************************
**************************************************************************************
*** A: Lower bend point

*** Local macros
local fig bwrange		/*data set identifier*/
local bp lower 			/*bend point*/
local cov_controls "age_filing male black hearings musc mental"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs

* Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef=.
save "${output}/figures/`fig'_`bp'.dta", replace
restore

* Calculating the ratio of b to h from main analysis
*> Note: as we vary bandwidth (h), we use the ratio of bias bandwidth (b) to it from the main analysis
preserve
keep if imeabs<=1200 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls') 
local b_h = e(b_l) / e(h_l)
display `b_h'
restore

* Running the analysis in $10 increments
forvalues i=700(-10)130 {

	local j = `i' * `b_h'

	capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) h(`i') b(`j') kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

	local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
	local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

* Saving the regression output
regsave using "${output}/figures/`fig'_`bp'.dta", autoid detail(scalars) addvar(r_m, `b_rob', `se_rob') append
}

*** Creation of figure

* Cleaning the data for figure creation
use "${output}/figures/`fig'_`bp'.dta", clear
keep var coef stderr h_l
drop if var=="RD_Estimate"
drop if coef==.
compress

* Creation of confidence intervals
gen double coef_l = coef - 1.96*stderr
gen double coef_u = coef + 1.96*stderr

* Define bandwidth used and add restrictions
rename h_l bw

keep if bw>=130
keep if bw<=700

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Bandwidth"
local yaxis "Estimate"
local xlabel "200(200)600"
local ylabel "-1(0.5)0.2"
local main_estimate 195
local labelsizes "medlarge"

#delimit ;
twoway (line coef bw, lc(black) lp(solid) lw(thick))
 (line coef_l bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 (line coef_u bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 , yline(0) xline(`main_estimate', lc(black) lp(longdash) lw(medium)) scheme(s1mono) 
 xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
 ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
 text(0.2 200 "{&larr}", placement(e) size(large)) 
 text(0.2 235 "Bandwidth used" "for main estimate", placement(e) justification(left) size(medium)) 
 title(`paneltitle', linegap(2) margin(b=4)) 
 legend(off) ;
#delimit cr
graph export "${output}/figures/Fig3A.eps", replace

**************************************************************************************
**************************************************************************************
*** B: Family bend point

*** Local macros
local fig bwrange		/*data set identifier*/
local bp fm 			/*bend point*/
local cov_controls "age_filing male black hearings musc mental"

use "${dta}/`bp'_main_sample.dta", clear

*** Creation of figure inputs

* Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef=.
save "${output}/figures/`fig'_`bp'.dta", replace
restore

* Calculating the ratio of b to h from main analysis
*> Note: as we vary bandwidth (h), we use the ratio of bias bandwidth (b) to it from the main analysis
preserve
keep if imeabs<=1200 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls') 
local b_h = e(b_l) / e(h_l)
display `b_h'
restore

* Running the analysis in $10 increments
forvalues i=1200(-10)100 {

	local j = `i' * `b_h'

	capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) h(`i') b(`j') kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

	local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
	local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

* Saving the regression output
regsave using "${output}/figures/`fig'_`bp'.dta", autoid detail(scalars) addvar(r_m, `b_rob', `se_rob') append
}

*** Creation of figure

* Cleaning the data for figure creation
use "${output}/figures/`fig'_`bp'.dta", clear
keep var coef stderr h_l
drop if var=="RD_Estimate"
drop if coef==.
compress

* Creation of confidence intervals
gen double coef_l = coef - 1.96*stderr
gen double coef_u = coef + 1.96*stderr

* Define bandwidth used and add restrictions
rename h_l bw

keep if bw>=330
keep if bw<=1200

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Bandwidth"
local yaxis "Estimate"
local xlabel "400(200)1200"
local ylabel "-0.8(0.4)0.2"
local main_estimate 600
local labelsizes medlarge

#delimit ;
twoway (line coef bw, lc(black) lp(solid) lw(thick))
 (line coef_l bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 (line coef_u bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 , yline(0) xline(`main_estimate', lc(black) lp(longdash) lw(medium)) scheme(s1mono) 
 xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
 ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
 text(0.1 605 "{&larr}", placement(e) size(large)) 
 text(0.1 645 "Bandwidth used" "for main estimate", placement(e) justification(left) size(medium))  
 title(`paneltitle', linegap(2) margin(b=4)) 
 legend(off) ;
#delimit cr
graph export "${output}/figures/Fig3B.eps", replace

**************************************************************************************
**************************************************************************************
*** C: Upper bend point

*** Local macros
local fig bwrange		/*data set identifier*/
local bp upper			/*bend point*/
local cov_controls "age_filing male black hearings musc mental"

use "${dta}/`bp'_main_sample.dta", clear

*** Creation of figure inputs

* Creating a data set to store the regression results 
preserve
clear
set obs 1
gen coef=.
save "${output}/figures/`fig'_`bp'.dta", replace
restore

* Calculating the ratio of b to h from main analysis
*> Note: as we vary bandwidth (h), we use the ratio of bias bandwidth (b) to it from the main analysis
preserve
keep if imeabs<=1200 
capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls') 
local b_h = e(b_l) / e(h_l)
display `b_h'
restore

* Running the analysis in $10 increments
forvalues i=1200(-10)430 {

	local j = `i' * `b_h'

	capture noisily rdrobust d_year14 ime, c(0) deriv(1) fuzzy(fmax) p(1) q(2) h(`i') b(`j') kernel(uni) scaleregul(0) vce(nn 3) covs(`cov_controls')

	local b_rob = e(tau_bc) * (1/12) * 1000  * 100 
	local se_rob = e(se_tau_rb) * (1/12) * 1000  * 100 

* Saving the regression output
regsave using "${output}/figures/`fig'_`bp'.dta", autoid detail(scalars) addvar(r_m, `b_rob', `se_rob') append
}

*** Creation of figure

* Cleaning the data for figure creation
use "${output}/figures/`fig'_`bp'.dta", clear
keep var coef stderr h_l
drop if var=="RD_Estimate"
drop if coef==.
compress

* Creation of confidence intervals
gen double coef_l = coef - 1.96*stderr
gen double coef_u = coef + 1.96*stderr

* Define bandwidth used and add restrictions
rename h_l bw

keep if bw>=430
keep if bw<=1200

* Figure settings
local paneltitle "C: Upper bend point"
local xaxis "Bandwidth"
local yaxis "Estimate"
local xlabel "600(200)1200"
local ylabel ""
local main_estimate 589
local labelsizes medlarge

#delimit ;
twoway (line coef bw, lc(black) lp(solid) lw(thick))
 (line coef_l bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 (line coef_u bw, color(gs12%50) lp(shortdash) lc(black) lw(medthick))
 , yline(0) xline(`main_estimate', lc(black) lp(longdash) lw(medium)) scheme(s1mono) 
 xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
 ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
 text(0.95 594 "{&larr}", placement(e) size(large))
 text(0.95 636 "Bandwidth used" "for main estimate", placement(e) justification(left) size(medium)) 
 title(`paneltitle', linegap(2) margin(b=4)) 
 legend(off) ;
#delimit cr
graph export "${output}/figures/Fig3C.eps", replace


**************************************************************************************
**************************************************************************************
******************************** APPENDIX FIGURES ************************************
**************************************************************************************
**************************************************************************************

**************************************************************************************
************************************ Figure B2 ***************************************
**************************************************************************************
*** Distribution of covariates around the lower bend point. (Produces Panels A-F.)

*** Local macros
local sample main
local bp lower 		/*bend point*/
local fig covar		/*data set identifier*/
local restrict 500	/*symmetric range for assignment variable*/

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
cap drop dens
gen dens=1

local dep_set1 = "age_filing"
foreach dep of local dep_set1 {

* Starting with the first covariate (age_filing) to create data set
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

* Add the other variables by merging
local dep_set2 = "male black hearings musc mental"
foreach dep of local dep_set2 {
		
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
merge 1:1 rdplot_id using "${output}/figures/`fig'_`bp'_`restrict'.dta"
sort rdplot_id
drop _merge
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

**************************************************************************************
*** Creation of figures
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

*****************
*** A: Age at filing
local dep "age_filing"
local paneltitle "A: Average Age when Starting DI"
local xaxis "Distance from Bend Point ($)"
local yaxis "Average age"
local xlabel "-400(200)400"
local ylabel "46(1)49"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.0fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2A.eps", replace

*****************
*** B: Males
local dep "male"
local paneltitle "B: Fraction of Males"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-400(200)400"
local ylabel "0.15(0.05)0.30"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2B.eps", replace

*****************
*** C: Black
local dep "black"
local paneltitle "C: Fraction Black"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-400(200)400"
local ylabel "0.09(0.02)0.15"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2C.eps", replace

*****************
*** D: Hearings
local dep "hearings"
local paneltitle "D: Fraction of Hearings Allowance"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-400(200)400"
local ylabel "0.24(0.03)0.345"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2D.eps", replace

*****************
*** E: Mental
local dep "mental"
local paneltitle "E: Fraction with Mental Disorders"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-400(200)400"
local ylabel "0.20(0.02)0.26"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2E.eps", replace

*****************
*** F: Musculo
local dep "musc"
local paneltitle "F: Fraction with Musculoskeletal Conditions"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-400(200)400"
local ylabel "0.28(0.02)0.33"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB2F.eps", replace


**************************************************************************************
************************************ Figure B3 ***************************************
**************************************************************************************
*** Distribution of covariates around the family bend point. (Produces Panels A-F.)

*** Local macros
local sample main
local bp fm 		/*bend point*/
local fig covar		/*data set identifier*/
local restrict 700	/*symmetric range for assignment variable*/

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
cap drop dens
gen dens=1

local dep_set1 = "age_filing"
foreach dep of local dep_set1 {

* Starting with the first covariate (age_filing) to create data set
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

* Add the other variables by merging
local dep_set2 = "male black hearings musc mental"
foreach dep of local dep_set2 {
		
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
merge 1:1 rdplot_id using "${output}/figures/`fig'_`bp'_`restrict'.dta"
sort rdplot_id
drop _merge
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

**************************************************************************************
*** Creation of figures
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

*****************
*** A: Age at filing
local dep "age_filing"
local paneltitle "A: Average Age when Starting DI"
local xaxis "Distance from Bend Point ($)"
local yaxis "Average age"
local xlabel "-600(200)600"
local ylabel "40(1)42"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.0fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3A.eps", replace

*****************
*** B: Males
local dep "male"
local paneltitle "B: Fraction of Males"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.40(0.10)0.60"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3B.eps", replace

*****************
*** C: Black
local dep "black"
local paneltitle "C: Fraction Black"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.14(0.02)0.19"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3C.eps", replace

*****************
*** D: Hearings
local dep "hearings"
local paneltitle "D: Fraction of Hearings Allowance"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.20(0.02)0.24"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3D.eps", replace

*****************
*** E: Mental
local dep "mental"
local paneltitle "E: Fraction with Mental Disorders"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.26(0.04)0.34"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3E.eps", replace

*****************
*** F: Musculo
local dep "musc"
local paneltitle "F: Fraction with Musculoskeletal Conditions"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.20(0.02)0.24"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB3F.eps", replace


**************************************************************************************
************************************ Figure B4 ***************************************
**************************************************************************************
*** Distribution of covariates around the upper bend point. (Produces Panels A-F.)

*** Local macros
local sample main
local bp upper		/*bend point*/
local fig covar		/*data set identifier*/
local restrict 700	/*symmetric range for assignment variable*/

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
cap drop dens
gen dens=1

local dep_set1 = "age_filing"
foreach dep of local dep_set1 {

* Starting with the first covariate (age_filing) to create data set
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

* Add the other variables by merging
local dep_set2 = "male black hearings musc mental"
foreach dep of local dep_set2 {
		
preserve 
keep if imeabs<=`restrict' 
capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
merge 1:1 rdplot_id using "${output}/figures/`fig'_`bp'_`restrict'.dta"
sort rdplot_id
drop _merge
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore
}

**************************************************************************************
*** Creation of figures
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

*****************
*** A: Age at filing
local dep "age_filing"
local paneltitle "A: Average Age when Starting DI"
local xaxis "Distance from Bend Point ($)"
local yaxis "Average age"
local xlabel "-600(200)600"
local ylabel "50(1)52"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.0fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4A.eps", replace

*****************
*** B: Males
local dep "male"
local paneltitle "B: Fraction of Males"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.68(0.02)0.76"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4B.eps", replace

*****************
*** C: Black
local dep "black"
local paneltitle "C: Fraction Black"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.11(0.01)0.13"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4C.eps", replace

*****************
*** D: Hearings
local dep "hearings"
local paneltitle "D: Fraction of Hearings Allowance"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.23(0.02)0.27"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4D.eps", replace

*****************
*** E: Mental
local dep "mental"
local paneltitle "E: Fraction with Mental Disorders"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.16(0.01)0.18"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4E.eps", replace

*****************
*** F: Musculo
local dep "musc"
local paneltitle "F: Fraction with Musculoskeletal Conditions"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction of sample"
local xlabel "-600(200)600"
local ylabel "0.28(0.01)0.305"
local labelsizes large

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(navy)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(cranberry) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB4F.eps", replace


**************************************************************************************
************************************ Figure B5 ***************************************
**************************************************************************************
*** Fraction of DI beneficiaries who received SSI [removed from sample]
*** (Produces Panels A, B & C.)

local sample ssi

**************************************************************************************
*** A: Lower bend point

*** Local macros
local fig ssi_dens	/*data set identifier*/
local bp lower		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "ssi"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction receiving SSI"
local xlabel "-400(200)400"
local ylabel "0.5(0.1)0.8"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB5A.eps", replace

**************************************************************************************
*** B: Family bend point

*** Local macros
local fig ssi_dens	/*data set identifier*/
local bp fm			/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "ssi"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction receiving SSI"
local xlabel "-600(200)600"
local ylabel "0.35(0.05)0.55"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB5B.eps", replace

**************************************************************************************
*** C: Upper bend point

*** Local macros
local fig ssi_dens	/*data set identifier*/
local bp upper	/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "ssi"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Upper bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Fraction receiving SSI"
local xlabel "-600(200)600"
local ylabel ""
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB5C.eps", replace


**************************************************************************************
************************************ Figure B6 ***************************************
**************************************************************************************
*** Actual DI payments around the bend points. (Produces Panels A, B & C.)

local sample main

**************************************************************************************
*** A: Lower bend point

*** Local macros
local fig payments	/*data set identifier*/
local bp lower 		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "fmax"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "DI payments ($)"
local xlabel "-400(200)400"
local ylabel "300(200)900"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB6A.eps", replace


**************************************************************************************
*** B: Family bend point

*** Local macros
local fig payments	/*data set identifier*/
local bp fm 		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "fmax"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "DI payments ($)"
local xlabel "-600(200)600"
local ylabel "1300(200)2100"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB6B.eps", replace


**************************************************************************************
*** C: Upper bend point

*** Local macros
local fig payments	/*data set identifier*/
local bp upper		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "fmax"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Upper bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "DI payments ($)"
local xlabel "-600(200)600"
local ylabel "2000(100)2300"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB6C.eps", replace


**************************************************************************************
************************************ Figure B7 ***************************************
**************************************************************************************
*** Placebo mortality rates for beneficiaries without dependents around the family
*** bend point

*** Local macros
local sample nondepend
local fig nondep	/*data set identifier*/
local bp fm			/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Sharp RK with covariates - no first stage in DI payments

use "${dta}/fm_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "3.7(0.3)4.65"
local labelsizes large

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB7.eps", replace


**************************************************************************************
************************************ Figure B8 ***************************************
**************************************************************************************
*** Predicted mortality rates for beneficiaries based on available covariates
*** (Produces Panels A, B & C.)

local sample main

**************************************************************************************
*** A: Lower bend point

*** Local macros
local fig predict	/*data set identifier*/
local bp lower 		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "covind_full"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Create covariate index
xtile age_q = age_filing, nq(10)
reg d_year14 i.age_q male black hearings musc mental
predict covind_full

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Predicted mortality rate (%)"
local xlabel "-400(200)400"
local ylabel "3.4(0.3)4.4"
local labelsizes large

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB8A.eps", replace


**************************************************************************************
*** B: Family bend point

*** Local macros
local fig predict	/*data set identifier*/
local bp fm 		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "covind_full"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Create covariate index
xtile age_q = age_filing, nq(10)
reg d_year14 i.age_q male black hearings musc mental
predict covind_full

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Predicted mortality rate (%)"
local xlabel "-600(200)600"
local ylabel ""
local labelsizes large

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB8B.eps", replace


**************************************************************************************
*** C: Upper bend point

*** Local macros
local fig predict	/*data set identifier*/
local bp upper 		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "covind_full"

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Create covariate index
xtile age_q = age_filing, nq(10)
reg d_year14 i.age_q male black hearings musc mental
predict covind_full

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace
restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Upper bend point"
local yaxis "Mortality rate (%)"
local xaxis "Distance from Bend Point ($)"
local xlabel "-600(200)600"
local ylabel "4.4(0.2)5.1"
local labelsizes large

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB8C.eps", replace


**************************************************************************************
************************************ Figure B9 ***************************************
**************************************************************************************
*** Placebo mortality rates for non-beneficiaries. (Produces Panels A, B & C.)

local sample nonbenef

**************************************************************************************
*** A: Lower bend point

*** Local macros
local fig nonbenef	/*data set identifier*/
local bp lower		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "mort"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-400(200)400"
local ylabel "0.20(0.05)0.35"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB9A.eps", replace


**************************************************************************************
*** B: Family bend point

*** Local macros
local fig nonbenef	/*data set identifier*/
local bp fm			/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "mort"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "0.20(0.05)0.35"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB9B.eps", replace


*************************************************************************************
*** C: Upper bend point

*** Local macros
local fig nonbenef	/*data set identifier*/
local bp upper		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "mort"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Upper bend point"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "0.25(0.05)0.40"
local labelsizes medium

#delimit ;
twoway	(scatter `dep'_mean `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.2fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB9C.eps", replace


**************************************************************************************
************************************ Figure B10 **************************************
**************************************************************************************
*** Mortality estimates for placebo kink locations. (Produces Panels A, B & C.)

local sample main

**************************************************************************************
*** A: Lower bend point

*** Local macros
local bp lower
tempvar v 
local h_bw = 194.480
local b_bw = 343.788
local max_bw = max(`h_bw',`b_bw')

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
gen `v' = ime

keep d_year14 `v'

local upper_bound = 1200
gen g = runiform(`max_bw', `upper_bound'-`max_bw')

capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) h(`h_bw') b(`b_bw') vce(nn 3)

local actual_t_cl = e(tau_cl)/e(se_tau_cl)
local actual_t_rb = e(tau_bc)/e(se_tau_rb)

gen pt_t_cl=0
gen pt_t_rb=0
gen pt_h=0
gen pt_b=0

gen shift_kink = 0

gen v_orig = `v'

forvalues i=1/200 {
	local gr=g[`i']
	qui replace `v' = v_orig - `gr'	
	qui replace shift_kink = `gr' if _n==`i'	

	** Placebo regressions **
	capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) vce(nn 3)
	
	qui replace	pt_t_cl = e(tau_cl)/e(se_tau_cl) if _n==`i'
	qui replace pt_t_rb = e(tau_bc)/e(se_tau_rb) if _n==`i'
	
	qui replace pt_h = e(h_l)
	qui replace pt_b = e(b_l)
	
				}

keep if pt_t_cl != 0
keep pt_t_cl pt_t_rb shift_kink pt_h pt_b

save "${output}/figures/pt_`bp'_random_bw.dta", replace

di "Actual conventional t-stat is `actual_t_cl'" 
di "Actual robust t-stat is `actual_t_rb'" 

*** Creation of figure
use "${output}/figures/pt_`bp'_random_bw.dta", replace

cumul pt_t_rb, gen(cum_t_rb)
sort pt_t_rb

* Figure settings
local paneltitle "A: Lower bend point"
local xaxis "Placebo t-statistics"
local yaxis "Empirical CDF"
local xlabel "-2(1)3"
local ylabel ""
local labelsizes medium

local true "2.62"

#delimit;
twoway (line cum_t_rb pt_t_rb, sort lpattern(solid) lcolor(black) lw(medthick)), 
	scheme(s1mono) xline(`true', lcolor(cranberry) lwidth(medium)) 
	xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1))
	ytitle(`yaxis', size(`labelsizes'))	ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))	
	text(0.8 2.35 "{&rarr}", placement(e) size(large)) 
	text(0.8 1.65 "Actual" "t-statistic", placement(e) justification(right) size(medium))
	title(`paneltitle', linegap(2) margin(b=4)) 
	legend(off) name(`bp'_`outcome', replace);
#delimit cr
graph export "${output}/figures/FigB10A.eps", replace


**************************************************************************************
*** B: Family bend point

*** Local macros
local bp fm
tempvar v 
local h_bw = 599.573
local b_bw = 1199.965
local max_bw = max(`h_bw',`b_bw')

use "${dta}/`bp'_main_sample.dta", clear

*** Creation of figure inputs
gen `v' = ime

keep d_year14 `v'

local upper_bound = 3500
gen g = runiform(`max_bw', `upper_bound'-`max_bw')

capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) h(`h_bw') b(`b_bw') vce(nn 3)

local actual_t_cl = e(tau_cl)/e(se_tau_cl)
local actual_t_rb = e(tau_bc)/e(se_tau_rb)

gen pt_t_cl=0
gen pt_t_rb=0
gen pt_h=0
gen pt_b=0

gen shift_kink = 0

gen v_orig = `v'

forvalues i=1/200 {
	local gr=g[`i']
	qui replace `v' = v_orig - `gr'	
	qui replace shift_kink = `gr' if _n==`i'	

	** Placebo regressions **
	capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) vce(nn 3)
	
	qui replace	pt_t_cl = e(tau_cl)/e(se_tau_cl) if _n==`i'
	qui replace pt_t_rb = e(tau_bc)/e(se_tau_rb) if _n==`i'
	
	qui replace pt_h = e(h_l)
	qui replace pt_b = e(b_l)
	
				}

keep if pt_t_cl != 0
keep pt_t_cl pt_t_rb shift_kink pt_h pt_b

save "${output}/figures/pt_`bp'_random_bw.dta", replace

di "Actual conventional t-stat is `actual_t_cl'" 
di "Actual robust t-stat is `actual_t_rb'" 

*** Creation of figure
use "${output}/figures/pt_`bp'_random_bw.dta", replace

cumul pt_t_rb, gen(cum_t_rb)
sort pt_t_rb


* Figure settings
local paneltitle "C: Upper bend point"
local xaxis "Placebo t-statistics"
local yaxis "Empirical CDF"
local xlabel "-2(1)1.5""
local ylabel ""
local labelsizes medium

local true "3.21"

#delimit;
twoway (line cum_t_rb pt_t_rb, sort lpattern(solid) lcolor(black) lw(medthick)), 
	scheme(s1mono) xline(`true', lcolor(cranberry) lwidth(medium)) 
	xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1))
	ytitle(`yaxis', size(`labelsizes'))	ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))	
	text(0.8 2.9 "{&rarr}", placement(e) size(large)) 
	text(0.8 1.95 "Actual" "t-statistic", placement(e) justification(right) size(medium))
	title(`paneltitle', linegap(2) margin(b=4)) 
	legend(off) name(`bp'_`outcome', replace);
#delimit cr
graph export "${output}/figures/FigB10B.eps", replace

#delimit;
twoway (line cum_t_rb pt_t_rb, sort lpattern(solid) lcolor(black) lw(medthick)), 
	scheme(s1mono) xline(`true', lcolor(cranberry) lwidth(medium)) 
	xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1))
	ytitle(`yaxis', size(`labelsizes'))	ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))	

	title(`paneltitle', linegap(2) margin(b=4)) 
	legend(off) name(`bp'_`outcome', replace);
#delimit cr	
graph export "${output}/figures/FigB10B.eps", replace


**************************************************************************************
*** C: Upper bend point

*** Local macros
local bp upper
tempvar v 
local h_bw = 588.943
local b_bw = 1141.986
local max_bw = max(`h_bw',`b_bw')

use "${dta}/`bp'_`sample'_sample.dta", clear

*** Creation of figure inputs
gen `v' = ime

keep d_year14 `v'

local upper_bound = 3000
gen g = runiform(`max_bw', `upper_bound'-`max_bw')

capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) h(`h_bw') b(`b_bw') vce(nn 3)

local actual_t_cl = e(tau_cl)/e(se_tau_cl)
local actual_t_rb = e(tau_bc)/e(se_tau_rb)

gen pt_t_cl=0
gen pt_t_rb=0
gen pt_h=0
gen pt_b=0

gen shift_kink = 0

gen v_orig = `v'

forvalues i=1/200 {
	local gr=g[`i']
	qui replace `v' = v_orig - `gr'	
	qui replace shift_kink = `gr' if _n==`i'	

	** Placebo regressions **
	capture noisily rdrobust d_year14 `v', c(0) deriv(1) p(1) q(2) kernel(uni) vce(nn 3)
	
	qui replace	pt_t_cl = e(tau_cl)/e(se_tau_cl) if _n==`i'
	qui replace pt_t_rb = e(tau_bc)/e(se_tau_rb) if _n==`i'
	
	qui replace pt_h = e(h_l)
	qui replace pt_b = e(b_l)
	
				}

keep if pt_t_cl != 0
keep pt_t_cl pt_t_rb shift_kink pt_h pt_b

save "${output}/figures/pt_`bp'_random_bw.dta", replace

di "Actual conventional t-stat is `actual_t_cl'" 
di "Actual robust t-stat is `actual_t_rb'" 

*** Creation of figure
use "${output}/figures/pt_`bp'_random_bw.dta", replace

cumul pt_t_rb, gen(cum_t_rb)
sort pt_t_rb

* Figure settings
local paneltitle "B: Family bend point"
local xaxis "Placebo t-statistics"
local yaxis "Empirical CDF"
local xlabel "-2(1)1.5"
local ylabel ""
local labelsizes medium

local true "0.81"

#delimit;
twoway (line cum_t_rb pt_t_rb, sort lpattern(solid) lcolor(black) lw(medthick)), 
	scheme(s1mono) xline(`true', lcolor(cranberry) lwidth(medium)) 
	xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1))
	ytitle(`yaxis', size(`labelsizes'))	ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))	
	text(0.95 0.63 "{&rarr}", placement(e) size(large)) 
	text(0.94 0.12 "Actual" "t-statistic", placement(e) justification(right) size(medium))
	title(`paneltitle', linegap(2) margin(b=4)) 
	legend(off) name(`bp'_`outcome', replace);
#delimit cr	
graph export "${output}/figures/FigB10C.eps", replace


**************************************************************************************
************************************ Figure B11 **************************************
**************************************************************************************
*** Mortality rates for DI/SSI beneficiaries after the waiting period. 
*** (Produces Panels A-F.)

local sample ssi

**************************************************************************************
*** A: Lower bend point - SSI only

*** Local macros
local fig ssi_mort	/*data set identifier*/
local bp lower 		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "A: Lower bend point - DI/SSI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-400(200)400"
local ylabel "2.0(0.5)3.0"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11A.eps", replace


**************************************************************************************
*** B: Lower bend point - DI + SSI

*** Local macros
local fig dissi_mort	/*data set identifier*/
local bp lower 		/*bend point*/
local restrict 500	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
*keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "B: Lower bend point - All DI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-400(200)400"
local ylabel "2.5(0.5)3.0"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11B.eps", replace


**************************************************************************************
*** C: Family bend point - SSI only

*** Local macros
local fig ssi_mort	/*data set identifier*/
local bp fm			/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "C: Family bend point - DI/SSI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "3.0(0.5)4.1"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11C.eps", replace


**************************************************************************************
*** D: Family bend point - DI + SSI

*** Local macros
local fig dissi_mort	/*data set identifier*/
local bp fm			/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
*keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "D: Family bend point - All DI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "3.0(0.5)4.0"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11C.eps", replace



**************************************************************************************
*** E: Upper bend point - SSI only

*** Local macros
local fig ssi_mort	/*data set identifier*/
local bp upper		/*bend point*/
local restrict 700	/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "E: Upper bend point - DI/SSI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "4.5(0.5)5.8"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11E.eps", replace


**************************************************************************************
*** F: Upper bend point - DI + SSI

*** Local macros
local fig dissi_mort	/*data set identifier*/
local bp upper			/*bend point*/
local restrict 700		/*symmetric range for assignment variable*/
local dep "d_year14"

*** Creation of figure inputs
use "${dta}/`bp'_`sample'_sample.dta", clear

* Sample restriction - only those receiving SSI
*keep if ssi==1

*** Creation of figure inputs
preserve 
keep if imeabs<=`restrict' 
cap drop dens
gen dens=1

capture rdplot `dep' ime, binselect(qs) genvars hide
capture collapse (mean) rdplot_mean_y rdplot_mean_bin (sum) dens, by(rdplot_id) fast
rename rdplot_mean_y `dep'_mean
rename rdplot_mean_bin `dep'_bin
rename dens `dep'_obs
save "${output}/figures/`fig'_`bp'_`restrict'.dta", replace

restore

*** Creation of figure
use "${output}/figures/`fig'_`bp'_`restrict'.dta", clear

* Figure settings
local paneltitle "F: Upper bend point - All DI beneficiaries"
local xaxis "Distance from Bend Point ($)"
local yaxis "Mortality rate (%)"
local xlabel "-600(200)600"
local ylabel "4.5(0.5)5.1"
local labelsizes medlarge

* Converting mortality rates to percentage points
gen `dep'_mean2 = `dep'_mean*100

#delimit ;
twoway	(scatter `dep'_mean2 `dep'_bin, msize(large) msymbol(circle) mcolor(black)), 
		scheme(s1mono) xline(0, lpattern(shortdash) lcolor(black) lwidth(thin))  
		xtitle(`xaxis', size(`labelsizes')) xlabel(`xlabel', labs(`labelsizes')) xsc(titlegap(1)) 
		ytitle(`yaxis', size(`labelsizes')) ylabel(`ylabel', labs(`labelsizes') angle(horizontal) format(%9.1fc)) ysc(titlegap(0))
		title(`paneltitle', linegap(2) margin(b=4)) 
		legend(off) name(`bp'_`dep', replace);
#delimit cr
graph export "${output}/figures/FigB11F.eps", replace

**************************************************************************************
log close